#include <stdio.h>
#include <cmath>
#include <algorithm>
#include <limits>
#include <float.h>

using namespace std;

extern "C" {
	int lsame_(char*, char*);
}

/* Table of constant values */

static int c__1 = 1;
static float c_b32 = 0.f;

double
pow_ri(float ap, int bp) {
	double pow, x;
	int n;
	unsigned long u;

	pow = 1;
	x = ap;
	n = bp;

	if (n != 0) {
		if (n < 0) {
			n = -n;
			x = 1 / x;
		}
		for (u = n;;) {
			if (u & 01)
				pow *= x;
			if (u >>= 1)
				x *= x;
			else
				break;
		}
	}
	return (pow);
}

double slamch_(char *cmach)
{
    /* System generated locals */
    int i__1;
    float ret_val;

    /* Local variables */
    static float t;
    int it;
    static float rnd, eps, base;
    int beta;
    static float emin, prec, emax;
    int imin, imax;
    int lrnd;
    static float rmin, rmax;
    float rmach;
    extern int lsame_(char *, char *);
    float small;
    static float sfmin;
    extern /* Subroutine */ int slamc2_(int *, int *, int *, float 
	    *, int *, float *, int *, float *);


/*  -- LAPACK auxiliary routine (version 3.1) -- */
/*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/*     November 2006 */

/*     .. Scalar Arguments .. */
/*     .. */

/*  Purpose */
/*  ======= */

/*  SLAMCH determines single precision machine parameters. */

/*  Arguments */
/*  ========= */

/*  CMACH   (input) CHARACTER*1 */
/*          Specifies the value to be returned by SLAMCH: */
/*          = 'E' or 'e',   SLAMCH := eps */
/*          = 'S' or 's ,   SLAMCH := sfmin */
/*          = 'B' or 'b',   SLAMCH := base */
/*          = 'P' or 'p',   SLAMCH := eps*base */
/*          = 'N' or 'n',   SLAMCH := t */
/*          = 'R' or 'r',   SLAMCH := rnd */
/*          = 'M' or 'm',   SLAMCH := emin */
/*          = 'U' or 'u',   SLAMCH := rmin */
/*          = 'L' or 'l',   SLAMCH := emax */
/*          = 'O' or 'o',   SLAMCH := rmax */

/*          where */

/*          eps   = relative machine precision */
/*          sfmin = safe minimum, such that 1/sfmin does not overflow */
/*          base  = base of the machine */
/*          prec  = eps*base */
/*          t     = number of (base) digits in the mantissa */
/*          rnd   = 1.0 when rounding occurs in addition, 0.0 otherwise */
/*          emin  = minimum exponent before (gradual) underflow */
/*          rmin  = underflow threshold - base**(emin-1) */
/*          emax  = largest exponent before overflow */
/*          rmax  = overflow threshold  - (base**emax)*(1-eps) */

/* ===================================================================== */

/*     .. Parameters .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */
/*     .. External Functions .. */
/*     .. */
/*     .. External Subroutines .. */
/*     .. */
/*     .. Save statement .. */
/*     .. */
/*     .. Data statements .. */
/*     .. */
/*     .. Executable Statements .. */
 
	
	base = numeric_limits<float>::radix;
	t = it = numeric_limits<float>::digits;
	lrnd = numeric_limits<float>::round_style == round_to_nearest;
	imin = numeric_limits<float>::min_exponent;
	imax = numeric_limits<float>::max_exponent;
	if (lrnd) {
	    rnd = 1.f;
	    eps = pow_ri(base, 1 - it) / 2;
	} else {
	    rnd = 0.f;
	    eps = pow_ri(base, 1 - it);
	}	
	rmax = numeric_limits<float>::max();
	prec = eps * base;
	emin = (float) imin;
	emax = (float) imax;
	sfmin = rmin;
	small = 1.f / rmax;
	if (small >= sfmin) {

/*           Use SMALL plus a bit, to avoid the possibility of rounding */
/*           causing overflow when computing  1/sfmin. */

	    sfmin = small * (eps + 1.f);
	}
    

    if (lsame_(cmach, "E")) {
	rmach = eps;
    } else if (lsame_(cmach, "S")) {
	rmach = sfmin;
    } else if (lsame_(cmach, "B")) {
	rmach = base;
    } else if (lsame_(cmach, "P")) {
	rmach = prec;
    } else if (lsame_(cmach, "N")) {
	rmach = t;
    } else if (lsame_(cmach, "R")) {
	rmach = rnd;
    } else if (lsame_(cmach, "M")) {
	rmach = emin;
    } else if (lsame_(cmach, "U")) {
	rmach = rmin;
    } else if (lsame_(cmach, "L")) {
	rmach = emax;
    } else if (lsame_(cmach, "O")) {
	rmach = rmax;
    }

    ret_val = rmach;
    
    return ret_val;

/*     End of SLAMCH */

} /* slamch_ */
